/**
 * HHNeuron.java
 *   Hodgkin-Huxley Neuron
 *
 * Robert McDougal
 * February 2007 - April 2007
 *
 * Changes:
 *  10 April 2007 - modified advance to store time history in double-valued time
 *                  and to interpolate lastSpikeTime
 */
import java.util.LinkedList ;
import java.util.Vector ;
import static java.lang.Math.exp ;
import static java.lang.Math.pow ;

// probably should redo synapse input part.  FIX.
public class HHNeuron extends SkeletalNeuronFramework {
    private double m , h , n ;
    private double theta = .5 ;
    private static final double SCALE_V = 100 ;
    private TheDiffEq diffEq = new TheDiffEq();
    private RungeKutta45 solver = null;
    private double absErrStep = .00001;
    private double minStep = .00001;
    private double FIRING_THRESHOLD = .5 ;
    private double maxH = .01 ;
    
    public HHNeuron( ) {
        String [ ] params = { "vNa" , "vK" , "vL" , "gNa" , "gK" , "gL" , "C" , "phi" ,
        "gsyn" ,"vSyn", "absErr/step" , "minStepSize" , "maxStepSize" } ;
        
        // initially at a random membrane potential between -1 and 1
        setPotential( 0. , 2 * Math.random( ) - 1 ) ;
        
        // m = minf , h = hinf , n = ninf
        m = .05 ;
        h = .6 ;
        n = .317 ;
        
        parameterNames = params ;
    }
    
    public String getName( ) {
        return "Hodgkin-Huxley Neuron" ;
    }
    
    public String getParam( String name ) {
        if ( name.equals( "vNa" ) ) {
            return "" + diffEq.vNa ;
        } else if ( name.equals( "vK" ) ) {
            return "" + diffEq.vK ;
        } else if ( name.equals( "vL" ) ) {
            return "" + diffEq.vL ;
        } else if ( name.equals( "gNa" ) ) {
            return ""  + diffEq.gNa ;
        } else if ( name.equals( "gK" ) ) {
            return "" + diffEq.gK ;
        } else if ( name.equals( "gL" ) ) {
            return "" + diffEq.gL ;
        } else if ( name.equals( "C" ) ) {
            return "" + diffEq.c ;
        } else if ( name.equals( "phi" ) ) {
            return "" + diffEq.phi ;
        } else if (name.equals( "vSyn") ) {
            return "" + diffEq.vSyn ;
        } else if ( name.equals( "absErr/step")) {
            return "" + absErrStep ;
        } else if ( name.equals( "minStepSize" ) ) {
            return "" + minStep ;
        } else if ( name.equals( "maxStepSize" ) ) {
            return "" + maxH ;
        } else if ( name.equals( "gsyn" ) ) {
            return "" + diffEq._gsyn ;
        } else {
            return null ;
        }
    }
    
    public void setParam( String name , String data ) {
        double value ;
        
        // all parameters are doubles
        try {
            value = Double.parseDouble( data ) ;
        } catch ( NumberFormatException e ) {
            // not a double; ignore the set request
            return ;
        }
        if ( name.equals( "vNa" ) ) {
            diffEq.vNa = value ;
        } else if ( name.equals( "vK" ) ) {
            diffEq.vK = value ;
        } else if ( name.equals( "vL" ) ) {
            diffEq.vL = value ;
        } else if ( name.equals( "gNa" ) ) {
            diffEq.gNa = value ;
        } else if ( name.equals( "gK" ) ) {
            diffEq.gK = value ;
        } else if ( name.equals( "gL" ) ) {
            diffEq.gL = value ;
        } else if ( name.equals( "C" ) ) {
            diffEq.c = value ;
        } else if ( name.equals( "phi" ) ) {
            diffEq.phi = value ;
        } else if (name.equals( "vSyn") ) {
            diffEq.vSyn = value ;
        } else if (name.equals("absErr/step")) {
            absErrStep = value ;
        } else if ( name.equals("minStepSize" ) ) {
            minStep = value ;
        } else if ( name.equals( "maxStepSize" ) ) {
            maxH = value ;
        } else if ( name.equals("gsyn" ) ) {
            diffEq._gsyn = value ;
        }
    }
    
    public INeuron clone( ) {
        HHNeuron n = new HHNeuron( ) ;
        n.setParam( getParams( ) ) ;
        return n ;
    }
    
    private class TheDiffEq implements IOdeRhs {
        public IAppliedCurrent iApp;
        
        private double alphaM( double v ) {
            if (v == 25) {
                return phi;
            }  else{
                return  phi * .1 * ( v - 25 ) / ( 1 - exp( - ( v - 25 ) * .1 ) );
            }
        }
        
        private double betaM( double v ) {
            return phi * 4 * exp( - v / 18. ) ;
        }
        
        private double alphaH( double v ) {
            return phi * .07 * exp( -v * .05 ) ;
        }
        
        private double betaH( double v ) {
            return phi * 1. / ( 1. + exp( - ( v - 30 ) * .1 ) ) ;
        }
        
        private double alphaN( double v ) {
            if (v == 10) {
                return phi/10;
            }   else {
                return phi * .01 * ( v - 10 ) / ( 1. - exp( - ( v - 10 ) * .1 ) ) ;
            }
        }
        
        private double betaN( double v ) {
            return phi * .125 * exp( - v / 80. ) ;
        }
        public double vNa = 115, vK = -12, vL = 11;
        public double gNa = 120, gK = 36, gL = .3;
        public double c = 1, phi = 1;
        public int time = 0;
        public double _gsyn = .25;
        public double vSyn = .01;
        
        public double[] eval(double t, double[] y) {
            double v = y[0];
            double m = y[1];
            double h = y[2];
            double n = y[3];
            double dv, dm, dn, dh;
            double temp;
            double temp2;
            double linPot, sj;
            
            
            dv = SCALE_V * iApp.amount(t) - gNa * h * (v-vNa) * pow(m, 3) - gK * (v - vK) * pow(n, 4) - gL * (v - vL);
            dm = alphaM( v ) * ( 1 - m ) - betaM( v ) * m;
            dh = alphaH( v ) * ( 1 - h ) - betaH( v ) * h;
            dn = alphaN( v ) * ( 1 - n ) - betaN( v ) * n;
            
            // add synaptic currents...
            
            temp = 0;
            for ( ISynapse s : syn ) {
                if (s.getPreCell().getPotential(t - s.getDelay()) > FIRING_THRESHOLD){
                    sj = 1;
                }   else {
                    sj = 0;
                }
                temp = temp + s.strength() * sj ;
            }
            
            dv += _gsyn * temp * (getPotential()-vSyn);
            
            
            
            
            double[] result = {dv, dm, dh, dn};
            return result;
        }
    }
    
    
    public void advance( ) {
        double vOld , tOld ;
        double v = getPotential() * SCALE_V;
        int i;
        double dt = .01;
        double[] y = null;
        
        // if never used before, init
        if (solver == null) {
            double[] y0 = {v , m, h, n};
            diffEq.iApp = Iapp;
            solver = new RungeKutta45(diffEq, y0, 0);
            solver.EPS_A = absErrStep;
            solver.MIN_STEP = minStep;
            solver.h = .001 ;
            solver.maxH = maxH ;
        }
        
//      The following is roughly equivalent to:
//        y = solver.advanceTo(time+1, Math.min(solver.h,.001));
//      except that it forces the computer to record potentials at each step
        
        diffEq.time = time;
        solver.h = Math.min( solver.h , .001 ) ;
        while ( solver.t0 < time + 1 ) {
            tOld = solver.t0 ;
            y = solver.step(solver.h);
            vOld = v ;
            v = y[0] / SCALE_V ;
            
            // interpolate to find time of last spike
            if ( vOld < FIRING_THRESHOLD && v >= FIRING_THRESHOLD ) {
                double frac = ( .5 - vOld ) / ( v - vOld ) ;
                lastSpike = tOld * ( 1 - frac ) + frac * solver.t0 ;
            }
            
            setPotential( solver.t0 , v ) ;
            if ( solver.h > solver.maxH ) {
                solver.h = solver.maxH ;
            }
            if ( solver.h < solver.MIN_STEP ) {
                solver.h = solver.MIN_STEP;
            }
        };
        
        time++;
    }
    
    
    // test driver
    // result looks blocky, but that's from only plotting 30 points...
    // but results are accurate since for each time step we calculate many smaller steps
    public static void main( String args [] ) {
        HHNeuron neuron = new HHNeuron( ) ;
        Curve c ;
        int i ;
        
        // set neuron's initial potential over threshhold
        neuron.setPotential( 0 , .21 ) ;
        
        // run for 30 time steps
        for ( i = 0 ; i < 30 ; i ++ ) {
            neuron.advance( ) ;
        }
        
        c = neuron.getCurve( ) ;
        c.isVisible = true ;
        
        Graph g = new Graph( ) ;
        g.addCurve( c ) ;
        g.setWindow( 0 , -.5 , 30 , 1.1 ) ;
        
        g.makeWindow( 640 , 480 , Graph.QUIT_ON_CLOSE , "Hodgkin-Huxley Neuron" ) ;
    }
}
